function track = trackerCart(MeshX,MeshY,NumSpace,NumAzimuth,NumPolar)
% function track = trackerCart(MeshX,MeshY,NumSpace,NumAzimuth,NumPolar)
%
% This function produces track length data for the MOC solution of a square
% region divided into rectilinear cells.
%
% References: 1) A. Hebert, "Applied Reactor Physics"
%             2) NumSpace. Smith, MIT NSE guest lecture notes
%             3) E. Sosnovsky, MOC pincell code
%
% J. Roberts,  03/14/2011
%
% Inputs:
%       MeshX       - Spatial divisions along x axis
%       MeshY       - Spatial division along y axis
%       NumSpace    - number of spatial divisions per azimuthal angle
%       NumAzimuth  - number of azimuthal angles (for 180 degrees)
%       NumPolar    - number of polar angles [1, 2, or 3]
%
% Output:
%        track       - Matlab datastructure with all tracking
%                      information as follows:
%        track.NumSurf       - number of surfaces surrounding geometry
%                              (unused for now)
%        track.NumReg        - number of regions
%        track.NumAzimuth    - number of azimuthal angles
%        track.NumSpace      - number of spatial spatial divisions
%        track.NumTracks     - number of trajectories
%        track.SurfPerim     - perimeter lengths of surfaces      [NumSurf]
%        track.Volume        - volumes of regions                  [NumReg]
%        track.SinPhi        - sines of tracking angles        [NumAzimuth]
%        track.CosPhi        - cosines of tracking angles      [NumAzimuth]
%        track.PertDist      - distance between track          [NumAzimuth]
%        track.IdxAngle      - angle index                      [NumTracks]
%        track.IdxSurfEnter  - index of entrance surface        [NumTracks]
%        track.IdxSurfExit   - index of exit surface            [NumTracks]
%        track.Enter         - track entrance coordinates    [NumTracks][2]
%        track.Exit          - track exit coordinates        [NumTracks][2]
%        track.Weight        - track weight                     [NumTracks]
%        track.NumSeg        - number of track segments         [NumTracks]
%        track.SegData().SegRegIdx() - regions crossed  [NumTracks][NumSeg]
%        track.SegData().SegLength() - segment length   [NumTracks][NumSeg]
%
% Note, regions can be crossed more than once (e.g. an annulus).  For each
% individual crossing, NumSeg is incremented. The structure as laid out
% follows Hebert and is rather general.
%
% Geometry assumed as follows:
%
%  ^ [7]...
%  | [4][5][6]
%  y [1][2][3]
%     x -> 

Pitch               = MeshX(end)-MeshX(1); % assuming square unit cell

track.NumSurf       = 4; % assumed for now
NumX                = length(MeshX)-1;
NumY                = length(MeshY)-1;
track.NumReg        = NumX*NumY;
track.NumPolar      = NumPolar;
track.NumAzimuth    = NumAzimuth;
track.NumSpace      = NumSpace;
track.SurfPerim     = Pitch*ones(track.NumSurf,1); % unused
track.NumSeg        = zeros(track.NumAzimuth,track.NumSpace);

% regions are ordered from inner-most to outer
track.Volume      = zeros(track.NumReg,1);


if (MeshX(1) ~= 0 || MeshY(1) ~=0)
    error('Keep it simple: set initial mesh boundaries at 0!')
end
if ( Pitch ~= (MeshY(end)-MeshY(1)) )
    error('We can only deal with SQUARE cells.  MeshX and MeshX are INCONSISTENT!')
end

% track.Volume(end) = Pitch^2 - pi*Radii(1)^2;
% for i = 1:track.NumReg-1
%     track.Volume(i) = pi * ( Radii(i+1)^2 - Radii(i)^2 );
% end

%--------------------------------------------------------------------------
% COMPUTE AZIMUTHAL DIVISION AND SPATIAL INTERCEPTS
%--------------------------------------------------------------------------

if ( mod( NumAzimuth, 2 ) == 1 )
    NumAzimuth = NumAzimuth + 1;
    disp( ' Warning: We require NumAzimuth to be even. ' )
    disp([' Setting NumAzimuth = ',num2str(NumAzimuth) ] )
end

track.NumTracks = NumAzimuth * NumSpace;

deltaPhi = pi/NumAzimuth;

track.Enter = zeros( NumAzimuth, NumSpace, 2 );
track.Exit  = track.Enter;
track.IdxSurfEnter = zeros( NumAzimuth, NumSpace );
track.IdxSurfExitt = track.IdxSurfEnter;
phivec = zeros(NumAzimuth,1);

for m = 1:NumAzimuth
    % Preliminary azimuth
    if ( m <= NumAzimuth*0.5 )
        tempPhi = deltaPhi * (m-0.5);
    else
        tempPhi = deltaPhi * (NumAzimuth+0.5-m);
    end
    % Number of x and y intercepts
    %   easy approach: we have NumSpace to assign to the x or y axis.  For
    %   a given azimuth, the fraction given to x should be proportional to
    %   sin(x), since a steeper angle means more going through x.  Likewise
    %   for y axis and cos(x).  Hence, take fraction going through x to be
    %   sin(x)/(cos(x)+sin(x)) = tan(x)/( 1+tan(x)).  Fraction through y is
    %   just 1 minus the x fraction.
    tanPhi = tan(tempPhi);
    Nx = abs( ceil( NumSpace * tanPhi/(tanPhi+1.0) ) );
    Ny = NumSpace - Nx;
    % Actual Azimuth (either 1st or 2nd quadrant)
    %   Do to rounding, we have different actual phi's defined by the ratio
    %   of x to y intercepts.
    if ( m <= NumAzimuth*0.5 )
        phi = atan ( Nx / Ny );
    elseif (  m > NumAzimuth*0.5 ) % 2nd quadrant
        phi = pi - atan ( Nx / Ny );
    else
        error (' Error in computing actual azimuths... ');
    end

    phivec(m)=phi;
    
    % insert angle index and associated sine and cosine
    track.IdxAngle( m ) = m;
    track.SinPhi( m ) = sin(phi);
    track.CosPhi( m ) = cos(phi);
    
    % Horizontal and vertical steps
    dx = Pitch / Nx;
    dy = Pitch / Ny;

    % Perpendicular distance between tracks
    track.PertDist( m ) = abs(dx*sin(phi));%sqrt(dx^2+dy^2);%dy / sin(phi);
    
    % Calculating the track -to - corner distances
    x2cd = 0.5*dx;  % jar, why are these needed?
    y2cd = 0.5*dy;
    
    % Track entrance and exit locations and boundary IDs
    if ( 0.0 < phi && phi <= 0.5*pi ) % 1st quadrant
        % Entrance and exit coordinates
        track.Enter( m, 1:Ny, 1 )               = zeros( 1, Ny );                 % jar, all x intercepts are zero for tracks w/ nonzero y intercepts
        track.Enter( m, 1:Ny, 2 )               = (Pitch-y2cd ):(-dy):(0);            % jar, insert the y intercepts
        track.Enter( m, (Ny+1):NumSpace, 1 )    = (x2cd):(dx):(Pitch);                % insert the x intercepts for the rest
        track.Enter( m, (Ny+1):NumSpace, 2 )    = zeros( 1, Nx );                 % and give zero for  y int of non-zero x ints
        track.Exit( m, 1:Nx, 1 )                = (x2cd):(dx):(Pitch);                % and the same for exit...but now x first
        track.Exit( m, 1:Nx, 2 )                = Pitch*ones( 1, Nx );
        track.Exit( m, (Nx+1):NumSpace, 1 )     = Pitch*ones( 1, Ny );
        track.Exit( m, (Nx+1):NumSpace, 2 )     = (Pitch-y2cd):(-dy):(0);
        % Boundary IDs
        track.IdxSurfEnter( m, 1:Ny )           = 1*ones( 1, Ny );                % incident on left face [1]
        track.IdxSurfEnter( m, (Ny+1):NumSpace )= 3*ones( 1, Nx );                % on bottom [3]
        track.IdxSurfExit( m, 1:Nx )            = 4*ones( 1, Nx );                % exit from top [4]
        track.IdxSurfExit( m, (Nx+1):NumSpace ) = 2*ones( 1, Ny );                % from right [2]
    elseif ( 0.5*pi < phi && phi <= pi ) % 2nd quadrant
        % Entrance and exit coordinates
        track.Enter( m, 1:Nx, 1 )               = (x2cd):(dx):(Pitch);
        track.Enter( m, 1:Nx, 2 )               = zeros ( 1, Nx );
        track.Enter( m, (Nx+1):NumSpace, 1 )    = Pitch*ones( 1, Ny );
        track.Enter( m, (Nx+1):NumSpace, 2 )    = (y2cd):(dy):(Pitch);
        track.Exit( m, 1:Ny, 1 )                = zeros( 1, Ny );
        track.Exit( m, 1:Ny, 2 )                = (y2cd):(dy):(Pitch);
        track.Exit( m, (Ny+1):NumSpace, 1 )     = (x2cd):(dx):(Pitch);
        track.Exit( m, (Ny+1):NumSpace, 2 )     = Pitch*ones( 1, Nx );
        % Boundary IDs
        track.IdxSurfEnter( m, 1:Nx )           = 3*ones( 1, Nx );
        track.IdxSurfEnter( m, (Nx+1):NumSpace )= 2*ones( 1, Ny );
        track.IdxSurfExit( m, 1:Ny )            = 1*ones( 1, Ny );
        track.IdxSurfExit( m, (Ny+1):NumSpace ) = 4*ones( 1, Nx );
    else % Wrong direction
        error (' Error in assigning entrance/exit coordinates ... ');
    end
    
    %--------------------------------------------------------------------
    % Track segment lengths 
    %--------------------------------------------------------------------
    RegIdx = zeros((NumX+NumY-1),1); % initialize these
    Length = RegIdx;
    
    for k = 1:NumSpace
        
        % Track coordinates at entrance and exit
        Ix = track.Enter( m, k, 1 );
        Iy = track.Enter( m, k, 2 );
        Fx = track.Exit( m, k, 1 );
        Fy = track.Exit( m, k, 2 );
        
        % Calculating (total) track length (for verification/debug)
        TmpLength = sqrt( (Fx-Ix)^2 + (Fy-Iy)^2 );
        
        % Algorithm probably has to go like this:
        %  while ~exited
        %    have current location
        %    figure out if traveling along x (right for phi<pi/2, left for 
        %      phi>pi/2to next cell involves a
        %      delta y that would or would not intercept the horizontal
        %      mesh divider above me
        %    if I do intercept the mesh above me, then my new y is the y
        %      value of that mesh and dx is dy/tanP (tanP =dy/dx)
        %    else my new x value is the next vertical mesh and dy = dx*tanP
        %    in both cases, length = sqrt(dx^2 + dy^2) (or dy/sin)
        %

        
        TanPhi = track.SinPhi(m)/track.CosPhi(m);        
        % I'm starting at Ix, Iy.  What cell do I start in?
        
        if ( TanPhi>0)
            I = sum(Ix>=MeshX);
        else
            I = sum(Ix>MeshX);
        end
        J = sum(Iy>=MeshY);
            

        NotFinished = 1;
        TrackCount = 0;
        
%         if (m==6 && k==5)
%             lala = 1;
%         end
%         hold on
%         plot(Ix,Iy,'ro','MarkerSize',10)
        
        while (NotFinished)
            
            % Upon entering this loop, the indices I and J correspond to
            % the cell the track is about to enter.  This is direction
            % dependent in with respect to I (as we can go left or right)
            % but we always go up, so J is angle-independent.
            
            % Add a new track (initialized to zero)
            TrackCount = TrackCount + 1;

            % Distance to the next vertical mesh edge; depends on angle.
            % Use of TanPhi in this way assumes phi doesn't include pi/2.
            if (TanPhi>0)
                DistToNextX = MeshX(I+1)-Ix;
            else
                DistToNextX = Ix - MeshX(I);
            end
            % Distance to the next horizontal mesh edge.  Since we do only
            % 0 to pi, we are always pointed up, and we look only above.
            DistToNextY = MeshY(J+1)-Iy;
            
            % Record the region into which we're about to enter
            RegIdx(TrackCount) = I + (J-1)*NumX;
           
            % "temp" records the difference between the vertical distance
            % traveled at angle phi if we went all the way to the closest
            % vertical edge.  If this distance is greater than the distance
            % to the next closest horizontal edge, then we intercept that
            % horizontal edge.
            temp = abs(DistToNextX*TanPhi) - DistToNextY;
            if ( abs(temp)>1e-10 && temp>0 ) % I hit top;  good rule of thumb for eps?
                
                dy = DistToNextY; 
                Iy = MeshY(J+1); 
                J = J+1;
                dx = dy/TanPhi;  
                Ix = Ix + dx;
                Length(TrackCount) = dy/track.SinPhi(m);
                
            % If that distance is not greater, than we move to the next 
            % vertical edge.
            elseif ( abs(temp)>1e-10 && temp< 0) % I hit side
                
                dx = DistToNextX; 
                if (TanPhi>0)
                    Ix = MeshX(I+1); 
                    I = I+1;
                else
                    Ix = MeshX(I);
                    I = I-1;
                end                
                dy = dx*abs(TanPhi);   
                Iy = Iy + dy;
                Length(TrackCount) = dy/track.SinPhi(m);
                
            % If the distances are (very close to being) the same, then we
            % have intersected a corner.  The tolerance is 1e-14, which
            % worked for simple tests to enforce what should be corner
            % crossings.
            else
                
                %disp('Going through a corner!!')
                if (TanPhi>0)
                    Ix = MeshX(I+1); 
                    I = I+1;
                else
                    Ix = MeshX(I);
                    I = I-1;
                end 
                Iy = MeshY(J+1);
                J = J+1;              
                Length(TrackCount) = DistToNextY/track.SinPhi(m);
                
            end
            % plot(Ix,Iy,'rx', 'MarkerSize',10)

            
            % If my X or Y index corresponds to having left the system,
            % then I quit.
            if ( I == length(MeshX) || I==0 || J == length(MeshY) )
                NotFinished = 0;
                % Record the number of tracks, their indices and lengths
                track.NumSeg(m,k) = TrackCount;
                track.SegData(m,k).SegRegIdx = RegIdx(1:TrackCount);
                track.SegData(m,k).SegLength = Length(1:TrackCount);      
                % verify track length
                err = abs(TmpLength-sum(track.SegData(m,k).SegLength));
                if (err>1e-12)
                    error(['Lengths not matched: m=',num2str(m),...
                        ' k=',num2str(k), ' err=',num2str(err)])
                end
            end
             
        end

    end % iteration over spatial dimension
    
end % iteration over azimuths 


% Calculating azimuthal weights
track.Weight = zeros(NumAzimuth,1);
track.Weight(1) = (  pi+phivec(2)-phivec(NumAzimuth) ) / ( 2*pi );
track.Weight(NumAzimuth) = ...
    (  pi + phivec(1) - phivec(NumAzimuth-1) ) / ( 2*pi );
% All other m
for m = 2:(NumAzimuth-1)
    track.Weight(m) = ( phivec(m+1) - phivec(m-1) ) / ( 2*pi );
end
ddd=1;

V = zeros(track.NumReg,1);
k=1;
for j=1:NumY
    for i=1:NumX
        V(k)= (MeshX(i+1)-MeshX(i))*(MeshY(j+1)-MeshY(j));
        k=k+1;
    end
end
% Corrected segment lengths (volume conservation) do this latertrack.PertDist( m )
AppxVol = zeros(NumAzimuth,track.NumReg);
for m = 1:NumAzimuth
    for k = 1:NumSpace
        %TotalLength = sum(track.NumSeg(m,k));
        for z = 1:track.NumSeg(m,k)   
            AppxVol(m,track.SegData(m,k).SegRegIdx(z)) = ...
                AppxVol(m,track.SegData(m,k).SegRegIdx(z)) + ...
                track.PertDist( m ) * track.SegData(m,k).SegLength(z);        
        end
    end
end
for m = 1:NumAzimuth
    for k = 1:NumSpace
        for z = 1:track.NumSeg(m,k)   
            track.SegData(m,k).SegLength(z) = ...
                track.SegData(m,k).SegLength(z) * ...
                V(i) / AppxVol(m,track.SegData(m,k).SegRegIdx(z));        
        end
    end
end

%--------------------------------------------------------------------------
end % function tracker




